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ABSTRACT 

X  1 

We  compute  the  decay  rates  for  stationary  perturbations  of  Poiseuille  flow  in  channels 

/ 

and  pipes.  The  decay  rates  are  found  by  solving  eigenvalue  problems  of  ordinary  differential 
equations,  where  the  eigenvalues  give  the  rate  of  decay  for  the  perturbation.  A  two-point 
boundary  value  method  is  used  to  compute  the  eigenvalues  yielding  efficient  and  accurate 
calculations.  For  the  channel  flow  problem,  the  results  are  in  agreement  with  previous 
calculations, (e.g.  [3],  [4],  [5],  j7],  [12])  ffiowever,  the  problem  of  determining  the  rate  of 
decay  for  a  fluid  motion  in  a  pipe  has  not  been  considered  before.  We-prove  tha>  for  the 


3 


I 


Stokes  problem  in  a  pipe  the  eigenvalues,  governing  the  rate  of  decay,  are  complex.  Wec 
carry  out  computations  for  small  and  moderate  Reynolds  numbers,  also  high  Reynolds 
number  computations  were  done  to  show  the  effectiveness  of  this  method. 

AMS(MOS)  Subject  Classifications:  65L15,  76D07 

Key  Words:  Navier-Stokes,  eigenvalue  problem.  Poiseuille  flow,'  PASVA3,  Reynolds 
number,  asymptotic 
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COMPUTATION  OF  THE  EIGENVALUES  FOR  PERTURBATIONS  OF 
POISEUILLE  FLOW  USING  A  TWO-POINT  BOUNDARY  VALUE  METHOD 

Gerardo  A.  Ache' 

1.  Introduction 

In  this  paper  we  are  concerned  with  the  eigenvalue  problem  which  governs  the  rate 
of  decay  for  a  stationary  perturbation  of  Poiseuille  flow.  We  consider  two-dimensional 
viscous  motions  in  channels  and  axi-symmetric  viscous  flow  in  pipes.  We  assume  that 
the  difference  between  the  base  flow  and  Poiseuille  flow  decays  exponentially  downstream 
(or  upstream).  It  is  then  possible  to  seek  solutions  to  the  Navier-Stokes  equations,  far 
downstream  (or  upstream),  that  are  a  perturbation  to  the  Poiseuille  profile  and  that  decay 
exponentially  in  the  axial  direction.  The  equations  can  then  be  linearized  yielding  an 
ordinary  differential  eigenvalue  system  where  the  eigenvalues  determine  the  rate  of  decay 
for  the  stationary  perturbation. 

By  use  of  the  stream  function  formulation,  it  is  possible  to  reduce  the  eigenvalue  sys¬ 
tem  to  a  single  fourth-order  differential  eigenvalue  equation  for  the  decay  of  the  stationary 
perturbation.  In  the  two-dimensional  case  this  differential  equation  is  very  similar  to  the 
Orr-Sommerfeld  stability  equation.  Results  regarding  the  computation  of  these  eigen¬ 
values,  in  the  case  of  channel  flow,  have  been  presented  previously,  e.g.  [3], [4], [5], [7], [12], 
however  when  the  fluid  motion  is  considered  axially  symmetric,  the  problem  of  determining 
the  rate  of  decay  for  the  stationary  perturbation  has  not  been  considered  before. 
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To  compute  the  eigenvalues  several  methods  have  been  used,  for  example,  a  spectral 
method  by  Bramley  j3j,  and  Bramley  and  Dennis  [4  ,  an  initial  value  method  by  Bramley 
and  Dennis  [5],  a  singular  perturbation  method  by  Wilson  i  12],  etc.  All  of  these  works  have 
dealt  with  the  two-dimensional  channel  flow  problem,  and  with  the  exception  of  Bramley  [3] 
all  of  them  have  considered  the  Navier-Stokes  equations  in  the  stream  function  formulation. 
The  work  of  Bramley  and  Dennis  [4;  appears  to  be  the  most  complete  among  all  others 
mentioned  above.  They  used  an  extension  of  the  Orszag’s  spectral  method  presented  in 
[9].  In  [3]  Bramley  used  the  same  method  to  compute  the  eigenvalues  using  the  primitive 
formulation  for  the  Navier-Stokes  equations,  i.e.,  he  used  the  velocity  field  and  the  pressure 
instead  of  the  stream  function.  However,  the  results  of  Bramley  3’  are  not  in  complete 
agreement  with  [4] .  This  disagreement  may  be  the  result  of  using  the  wrong  number  of 
boundary  conditions.  There  are  two  disadvantages  reported  in  [4]  with  respect  to  spectral 
methods,  these  are  the  computation  of  spurious  eigenvalues  and  the  loss  of  accuracy  for 
high  Reynolds  numbers  computations,  especially  for  eigenvalues  with  negative  real  part. 
In  [5]  Bramley  and  Dennis  used  an  initial  value  method  to  compute  the  eigenvalues.  They 
said  that  using  this  method  it  is  possible  to  overcome  some  of  the  difficulties  of  spectral 
method.  However,  their  method  computes  only  eigenvalues,  and  not  eigenfunctions,  and 
only  one  at  a  time. 

In  this  paper  we  compute  the  eigenvalues  for  a  perturbation  of  Poiseuille  flow,  in  a 
channel  and  a  pipe,  using  a  more  accurate  and  efficient  method.  This  method  can  be 
described  in  two  steps,  first,  we  transform  the  eigenvalue  problem  into  an  equivalent  two- 


point  boundary  value  system,  then  we  numerically  solve  this  system  using  the  two-point 
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boundary  value  solver  DVCPR  from  the  IMSL  library.  For  the  two-dimensional  problem 
our  results  are  in  full  agreement  with  the  previous  work  of  Bramley  and  Dennis  [4],  for  small 
and  moderate  Reynolds  numbers,  and  with  the  asymptotic  prediction  of  Wilson  [12],  for 
eigenvalues  approaching  zero  at  high  Reynolds  numbers.  The  numerical  results  presented 
in  this  paper  show  the  efficiency  of  the  calculations  and  the  superiority  of  our  method  with 
respect  to  the  other  methods  mentioned  above.  These  eigenvalues  are  important  for  the 
derivation  of  boundary  conditions,  this  problem  was  considered  in  [l  . 


2.  Channel  Flow 

Given  a  semi-infinite  channel,  the  incompressible  Navier-Stokes  equations  are  given 


by 


du  du  dp  1  _2 

UTx+VTy  +  Tx-  RX  “  ’ 

(2-la) 

dv  dv  dp  1  _2 

u  —  +  v  -—  +  =  —  V  lv  . 

dx  dy  dy  R 

(2.16) 

du  dv 

ai  +  ai=°' 

(2.1c) 

with  boundary  conditions 

u  =  v  =  0  at  y  =  ±1  , 

(2*1*0 

u  =  F, (y)  ,  v  =  F2(y)  at  x  =  0  , 

(2.1c) 

where  F ,  is  a  profile  satisfying  F,  (-1)  =  0.  Finally  we  have  the  regularity  condition 

(u,r)  — ■  (u,0)  as  x  —  oc  .  (2.1/) 

Here  R  represents  the  Reynolds  number,  V2  =  d2  ,dx2  —  d2/dy2,  and  u  =  |(1  -  y2)  is  the 


Poiseuille  parabolic  profile. 


Since  we  are  only  interested  in  solutions  of  this  system  which  decay  exponentially,  far 


downstream  we  seek  an  asymptotic  solution  to  (2.1)  in  the  form 


u{x,y)  -  -(1  ~  S'2)  +  Wi  (y)  exp(-Ax)  , 


u(i,y)  =  Wr2(y)exp(-Ax)  , 


p(x,y)  =  -3R  lx~R  1q[y)  exp(-Ax)  -  C  , 


(2-2  a) 


(2.26) 


(2.2c) 


where  C  is  an  arbitrary  constant. 


Substituting  these  expressions  in  equations  (2.1a)  to  (2.1c),  and  neglecting  nonlinear 
terms,  we  have  that  (A,  W%, q)  satisfies  the  following  eigenvalue  system  of  ordinary 

differential  equations, 


__1  =  -/?{-(!  -  y*)\Wl  -  3 yW2}  -  A q  -  \H\\ 


=  AH7!  , 


(2.3a) 


(2.36) 


(2.3c) 


and  boundary  conditions 


W  i(±l)  =  H72(±l)  =  0  . 


(2.3d) 


We  are  interested  in  soKing  this  two-point  boundary  value  system  to  determine  the 
rate  of  decay  A.  For  x  large  (i.e..  downstream)  the  solution  to  (2.1)  may  be  represented  by 


2  _ ^ 

“  (x^y)  =  -  y2)  +  H'i,„(y)exp(-Anx)  , 

n 

ViX^y)  =  y"ir2.n(y)<“Xp(-Anx)  , 


(2.4a) 


(2.46) 
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and 


q{x,y)  =  -3R  lx  +  R  1  ^  gn(y)  exp(- Anx)  +  C  , 

n 

Re{ A„)  <  7?e(An+1)  . 


(2.4c) 


(2  Ad) 


The  sum  in  (2.4a)  to  (2.4c)  is  taken  over  all  the  eigenvalues  with  positive  real  part  (the 
eigenvalues  with  negative  real  part  can  be  used  to  construct  an  asymptotic  solution  to  (2.1) 
in  a  semi-infinite  channel  on  the  negative  x-axis)  and  (W'iifll  ^  2,n,  9,  An)  being  a  solution 
to  the  system  (2.3). 

We  may  eliminate  the  perturbed  pressure  q  from  the  system  (2.3)  by  differentiating 
(2.3a)  with  respect  to  y,  multiplying  equation  (2.3a)  by  —A,  adding  the  resulting  expres¬ 
sions  and  using  (2.36),  then  the  system  (2.3)  is  reduced  into  a  single  fourth-order  differen¬ 
tial  equation  in  W2  and  A  very  similar  to  the  Orr-Sommerfeld  equation  (for  simplicity  we 
replace  W2  by  W) 


d4W 

dyA 


,  d2W  .  ,  3 ,  7,,d2W 

2X2 -+-  A4W  -f  Ai2{-(1  -  y2)( 


dy‘ 


dyi 


A 2W)  +  3 W}  =  0  , 


(2.5  a) 


and  boundary  conditions 

AW 

W  =  — -  =0  at  y  =  ±1  .  (2.56) 

dy 

We  solved  this  equation  in  order  to  compare  our  method  of  solution  with  previous 
computations.  Equation  (2.5)  is  associated  with  the  stream  function  formulation  of  the 
Navier-Stokes  equations  (see  e.g.  4  ). 

When  the  Reynolds  number  R  is  equal  to  zero  we  may  find  an  explicit  solution  to  (2.5) 
and  therefore  to  (2.3).  This  solution  is  given  in  terms  of  the  Papkovitch-Fad'le  functions 
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for  the  biharmonic  equation  in  a  semi-infinite  strip  .  For  R  -  0  equation  (2.5)  becomes 


d4W  ,d2lV’  4 

+  2A2-t-^  +  A4W  -  0  , 


dy4 


dy- 


(2.6) 


and  boundary  conditions  (2.56).  This  equation  is  associated  with  the  biharmonic  equation 
when  solutions  of  the  form  w(x,y)  =  W{ y)exp(-Ax)  are  sought.  There  are  two  types  of 
solutions  for  (2.6),  one  even  and  the  other  odd,  which  are  given  by 


W  e(y)  =  (y  -  1)  sin(Ae(y  -f  1))  +  (y  1)  sin(Ae(y  -  1)) 


(2.7a) 


and 


W  °(y)  =  (y  -  1)  sin(A°(y  +  1))  -  (y  -*-  1)  sin(A°(y  -  1))  , 


(2.76) 


respectively.  To  satisfy  the  boundary  conditions  (2.56)  we  need  that  the  eigenvalues  Ae 
and  X°  satisfy  the  following  transcendental  equations, 

sin  (2Ae)  —  2Ae  =  0  ,  (2.8a) 

and 

sin  (2A°)  -  2AP  =  0  .  (2.86) 

respectively.  We  notice  that  if  A  is  an  eigenvalue  satisfying  (2.8)  so  is  -  A  and  also  ±A. 
also  the  function 

w{x,y)  =  y^qn»n(y)exp(-Anj)  .  (2-9) 

n 

where  the  sum  is  taken  over  the  eigenvalues  with  positive  real  part,  satisfies  the  biharmonic 
equation  in  a  semi-infinite  strip  with  transversal  boundary  conditions  w  =  du'  'dy  =  0  . 

In  section  4  we  discuss  a  method  to  solve  the  eigenvalue  problems  (2.3)  and  (2.5), 
where  the  solution  to  the  Stokes  problem  (2.G)  plays  an  important  role.  Also  in  section  5 
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we  present  numerical  results  involving  the  computation  of  these  eigenvalues  for  arbitrary 


Reynolds  numbers. 


3.  Axially  Symmetric  Flow 

In  this  section  we  consider  the  incompressible  Navier-Stokes  equations  in  cylindrical 
coordinates  and  dimensionless  form,  set  in  a  semi-infinite  pipe.  These  equations  can  be 


written  as, 


du  du  dp  1  .  u  \ 

u  — +  »-  +  -  =  -(V-U-  -3)  , 


dr 


dz  dr  R 


dv  dv  dp  1 

u  —  +  u—  +  —  =  —  V2u  , 
dr  dz  dz  R 


1  d  dv 

-,Tr(ru)+ai  =  a' 


and  boundary  conditions 


u  =  r  =  0  at  r  =  1  . 


(3.1a) 

(3.16) 

(3.1c) 


(3-lrf) 


Since  the  fluid  motion  is  symmetric  the  conditions  at  the  center  line  are 


dv 

u  =  —  =  0  at  r  =  0  , 
dr 


(3.1c) 


we  also  specify  the  entry  condition 


u  =  Fj(r)  ,v  =  F2(r)  at  2  =  0, 


(3.1/) 


finally  we  have  the  regularity  condition. 


(u,t>)  — *  (O.u)  as  2  — ►  00 
Here  V2  =  1/r  d/dr  (rd/dr  )  -  d2 /dz2  ,  and  v  =  1  -  r2. 


(3.1?) 


Similar  to  the  channel  flow  problem  we  seek  an  asymptotic  solution  to  (3.1)  in  the 


form, 


u(r,2 )  =  Wi[r)  exp(-A;r)  , 


(3.2a) 


v  (r,  z)  =  1  -  r2  -f  W'2(r)  exp(-  Az)  , 


(3.26) 


p  (r,  z)  =  —  4R  lz-R  ’g(r)  exp(- Az)  +  C  , 


(3.2c) 


with  C  an  arbitrary  constant. 

Substituting  these  expressions  in  (3.1a)  to  (3.1c),  and  neglecting  nonlinear  terms,  we 
obtain  an  eigenvalue  system  in  W\  ,  W2  .  q  and  A  similar  to  (2.3),  i.e.  , 


dr  r 


(3.3  a) 


=  -R{2rW1  +  A(l  -  r2)»'2}  -  -  ^  -  X2\V2  -  A q  , 
dr 2  f  dr 

^  =  R{1  -  r2)AH'i  -  A^  4-  A2Wj  , 
dr  dr 


(3.36) 


(3.3c) 


with  boundary  conditions, 


W  i(0)  =  — — — - (0)  =  0  and  H’,(l)  =  W'2(l)  =  0  . 
dr 


(3.3  d) 


By  solving  this  system  we  may  represent  the  solution  to  (2.1).  for  c  far  downstream, 
in  the  same  way  as  for  the  channel  flow  problem  i.e.,  by  expressions  similar  to  (2.4a)  to 


(2-4  d). 


Similar  to  the  channel  flow  problem  we  solved  the  system  (3.3)  numerically,  using  the 


solution  to  the  system  for  R  -  0  (the  Stokes  problem).  These  solutions  may  be  found  as 


follows.  When  R  =  0  the  system  (3.3)  becomes, 

=w3-?A, 

dr  r 

(3.4  a) 

W  -  \2W2  -  , 

dr 2  r  dr 

(3.46) 

dr  dr 

(3.4c) 

and  boundary  conditions  (3.3 d). 

To  find  an  explicit  solution  to  this  system  differentiate  (3.4c)  with  respect  to  r,  mul¬ 
tiply  (3.4a)  by  A2  and  (3.46)  by  A,  then  add  the  resulting  expressions  and  using  (3.4c)  we 
obtain  the  following  differential  equation  in  q, 


dll+1-f  +  x^  =  0.  (3.5) 

dr2  r  dr  v  ' 

This  equation  has  the  particular  solution, 

q[r)  =  J0{\r)  ,  (3.6) 


where  J o  is  the  Bessel  function  of  first  kind  and  order  zero.  This  is  the  only  solution,  up 
to  a  multiplicative  constant,  which  is  finite  when  r  is  equal  to  zero.  Substituting  q(r)  in 
(3.46)  we  obtain  the  following  differential  equation  in  IV2. 


d2W2  1  dWj2  x  2 


dr 2  t  dr 


-  A 2VV2  =  -AJo(Ar)  , 


with  boundary  conditions. 


dW  2 


dr 


(0)  =  W'2(l)  =  0  . 


This  two-point  boundary  value  problem  has  the  following  solution. 


W  2(r)  =  bJ0(*r)  -  -rJt(\r)  , 


(3.7a) 


(3.76) 


(3.8) 
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where  b  is  a  constant  to  be  determined  by  U;2(l)  =  0  .  Using  (3.4a)  then  IV  j  is  given  as, 


W  ,(r)  -  \[^rJ0(Xr)  -  (1  -  Afe^Ar)]  .  (3.9) 

Since  the  solutions  (3.8)  and  (3.9)  need  to  satisfy  the  boundary  conditions  (3.3d)  we  have 
that, 

AW,(l)  =  ^Jo(A)  -  (1  -A6)J,(A)  =0  ,  (3.10a) 

and 

W2(1)=6Jo(A)-^J1(A)=0.  (3.106) 

By  eliminating  b  from  these  two  equations  we  obtain, 

!^LEZj0(A)-U,(A)=0.  (3.11) 

2A  2 

To  solve  this  equation  we  observe  that  (3.11)  can  be  transformed,  after  some  algebraic 
manipulations,  into  the  following  non-linear  equation 

J  i(A)Jo(A)  -  ^(Jf( A)  +  J|( A))  =  0  .  (3.12) 


Lemma  3.1 

All  non-trivial  solutions  to  the  equation  (3.12),  that  is  those  with  A  =*  0  .  are  complex. 

Proof 

We  make  use  of  the  following  integral  formula  involving  the  Bessel  functions. 

J  xjf(x)dx  =  ^x2(J2(z)  -  J2(x))  -  zJq[x)Ji(t)  .  (3.13) 

10 
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a.  A.  A.  a.  A.  a.  A  .*  A.  .  A.  /.  /.  A  A  ’  a\A  A  AA  AAA  A  A  • 


This  formula  is  obtained  by  integrating  by  parts  Lommel’s  integral  formula  with  a  —  1 


(see  [2],  p  10).  Then  assuming  that  some  A  satisfying  (3.12)  is  real  we  get  the  contradiction 
that 

r *  i 

0  <  /  xJ,{x)2dx  =  -A2(J2(A)  -  J2(A))  -  AJ0(A)Ji(A)  =  0  , 

Jo  *• 

where  we  have  assumed,  without  loss  of  generality,  that  A  >  0  .  This  contradiction  proves 
the  lemma. 

It  is  known  that  the  Bessel  functions  have  the  following  asymptotic  representations, 
J  o(A)  =  ( — r)  1/2{cos(A  -~)-t  e-i(A)}  ,  (3.14a) 

7T  A  4 

and 

■,,{A)  =  (lt)''2{sin<A~  j>~£2<a»  •  (3-H4) 

where  cj(A),£:2(A)  can  be  regarded  as  small  order  terms  for  /?e(A)  large.  Substituting 
these  asymptotic  expressions  in  (3.12),  and  neglecting  £i(A)  and  52(A),  the  A's  can  be 
approximated  by  solutions  of  the  following  transcendental  equation, 

cos  (2A)  +  A  =  0  .  (3.15) 

Solutions  of  this  equation  can  be  used  as  initial  iterates  for  a  Newton  method  to  obtain 
the  roots  of  (3.12).  We  will  present  the  detail  of  the  computation  in  section  5. 

Similar  to  the  channel  flow  problem  we  used  the  solution  to  the  Stokes  problem  to 
numerically  solve  the  system  (3.3)  for  arbitrary  Reynolds  numbers. 

4.  Method  of  Solution 

In  this  section  we  discuss  a  method  to  solve  the  two-point  boundary  value  problems  for 
the  systems  (2.3),  (3.3)  and  equation  (2.5)  and  determine  their  eigenvalues.  The  method 
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can  be  described  as  follows,  first  the  eigenvalue  problem  is  transformed  into  an  equivalent 


two-point  boundary  value  first-order  system  of  the  form, 


Z'  =  A(A,  R,  t)Z  ,  a  <  t  <  b  , 


(4.1a) 


A'  =  0  , 


(4.16) 


with  boundary  conditions 


BaZ(a)  -  Bt.Z{6)  =  0 


(4.1c) 


Zfc(6)  =  1  , 


(4. Id) 


for  some  k  ,  1  <  k  <  N  .  Here  Z  =  (21,  Z2,  •  •  • ,  zn)  »  and  A,Ba,Bb  are  N  x  N  matrices. 


The  prime  in  (4.1a)  indicates  derivatives  with  respect  to  t.  The  condition  (4. Id)  has  to 


be  chosen  in  a  way  that  is  not  in  conflict  with  (4.1c).  In  the  case  of  systems  (2.3),  (3.3) 


and  equation  (2.5)  we  used  for  condition  (4. Id)  H’^l)  =  1.  W'{(1)  =  1,  and  IV"(1)  —  1. 


respectively.  To  numerically  solve  the  system  (4.1)  we  used  the  two-point  boundary  value 


system  code  DYCPR  from  the  IMSL  library  (this  code  is  also  known  as  PASVA3),  it  is 


a  standard  solver  for  first-order  system  of  ordinary  differential  equations  with  conditions 


at  two  end  points.  It  uses  a  variable  step,  with  an  automatic  criterion  to  select  a  non- 


uniform  grid,  and  a  variable  order  of  accuracy,  with  an  excellent  correspondence  between 


the  requested  tolerance  ( tol )  and  the  actual  global  error  in  the  numerical  solution,  (see 


[10).  Since  the  systems  considered  are  non-linear  we  used  as  initial  iterate  in  PASYA3  the 


solution  to  the  corresponding  Stokes  problem  {R  =  0)  ,  which  is  known  and  for  which  the 


boundary  conditions  are  satisfied,  to  generate  solutions  for  arbitrary  Reynolds  numbers. 


This  procedure  is  called  continuation.  The  disadvantages  of  this  method  are  that  we  have 


to  compute  the  eigenvalues  one  at  a  time,  however  the  eigenvalues  are  computed  together 
with  their  eigenfunctions.  Also  we  can  not  compute  the  solution  at  a  critical  Reynolds 
number,  i.e.  a  Reynolds  number  for  which  the  solution  to  the  system  is  not  isolated, 
however  it  is  possible  to  compute  solutions  for  Reynolds  numbers  close  to  a  critical  value. 
By  using  this  method  we  avoid  the  computation  of  spurious  eigenvalues  and  also  for  high 
Reynolds  numbers  the  solution  can  be  computed  in  a  very  accurate  manner. 

5.  Numerical  Results 

We  now  present  some  numerical  results  concerning  the  eigenvalue  problems  described 
in  §2  and  §3.  For  the  channel  flow  problem  our  results  are  in  full  agreement  with  the 
previous  computation  of  Bramley  and  Dennis  [4],  also  we  obtain  the  same  answer  by 
solving  the  system  (4.1)  corresponding  to  the  primitive  formulation,  i.e.,  equations  (2.3), 
and  the  system  (4.1)  corresponding  to  the  stream  function  formulation,  i.e..  equation  (2.5). 

We  present  about  the  same  number  of  results  as  in  '4].  We  present  our  results  graph¬ 
ically  concerning  the  behavior  of  these  eigenvalues  in  the  complex  plane.  We  regard  these 
eigenvalues  as  a  function  of  Reynolds  numbers  in  the  form  A(R)  =  a(R)  4-  id[R)  then  we 
plotted  the  pairs  (a(R), 3(R))  in  the  A  plane.  We  computed  even  and  odd  branches  of 
eigenvalues  and  eigenvalues  for  axi-symmetric  flow.  Note  that  for  the  axi-symmetric  flow 
there  is  no  associated  parity,  such  as  even  or  odd,  thus  when  we  refer  to  an  even  or  odd 
solution  it  means  an  even  or  odd  solution  to  the  channel  problem. 

For  the  channel  flow  problem  we  have  found  that  these  eigenvalues  behave  as  follows 
(see  figures  1.1  a  and  1.1  b).  Let  (A*(0),  Z*(l))  be  an  eigenpair  solution  for  R  —  0  where 
Re(Afc)  :<]  /?e(A*+1)  |  ,  k  =  1.2,---  .  For  0  <  R  <  R*  0  (A*(R).  Z*(<))  remains 
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complex,  its  conjugate  also  being  an  eigenpair  solution.  At  R  —  this  complex  solution 
and  its  conjugate  coalesce  on  the  real  axis  and  then,  for  R  larger  than  R* 0  they  split 
into  two  branches  of  real  solutions,  one  with  increasing  eigenvalues  and  the  other  with 
decreasing  eigenvalues.  For  solutions  with  positive  eigenvalues  the  increasing  real  solution 
(A*(/2),  Zk{t))  coalesces  with  the  decreasing  real  solution  (A*,.  j  (/?),  Zfc+i(f))  of  the  next 
branch  at  R  =  R*  j  then,  for  larger  values  of  R  these  two  solutions  split  into  a  complex 
solution  and  its  conjugate  and  again  they  remain  complex  for  R* ,  <  R  <  R*  2  ,  then 
at  R  =  R^2  the  complex  solution  and  its  conjugate  again  coalesce  on  the  real  axis  and 
the  cycle  is  repeated.  Notice  that  R^m,  m  =  1,2,--*  ,  are  the  critical  Reynolds  numbers. 
There  is  an  exception  to  this  cycling  phenomenon  the  first  even  branch  of  solution  remains 
complex  and  approaches  zero  as  R  increases.  Also  the  first  real  branch  of  even  and  odd 
decreasing  eigenvalues  approaches  zero  as  R  increases.  None  of  the  previous  researchers 
have  reported  the  second  type  of  critical  Reynolds  number  at  which  two  real  solutions 
coalesce  and  split  into  a  complex  solution  and  its  conjugate.  Only  Wilson  i  12]  pointed  out 
that  the  positive  increasing  real  eigenvalues  approach  a  fixed  value  on  the  real  axis.  As  can 
be  seen  from  Figures  1.1a  and  1.1b  the  cycles  of  the  eigenvalues  appear  to  be  converging 
to  a  value  on  the  real  axis. 

For  eigenvalues  with  negative  real  part  the  situation  is  different,  solutions  with  eigen¬ 
values  which  decrease  in  magnitude  either  approach  zero  as  R  increases,  but  in  a  much 
slower  rate  than  solutions  with  positive  real  eigenvalues,  or  they  approach  a  fixed  value. 
In  these  cases,  the  real  eigenvalues  increase  or  decrease  very  slowly  that  for  the  range  of 


WWW 


the  Reynolds  number  for  which  we  made  the  computations,  we  did  not  observe  the  coa¬ 
lescence  of  real  eigenvalues  with  negative  real  part. 


We  have  found  that  for  the  axi-symmetric  problem  the  solutions  behave  in  the  same 
manner  as  the  odd  solutions  for  the  channel  problem.  That  is,  the  solution  are  initially 
complex,  the  conjugate  pairs  coalesce  on  the  real  axis,  the  real  solutions,  with  positive 
eigenvalues,  coalesce,  becoming  complex,  etc. 

Figures  1.1a  and  1.16  show  the  eigenvalues  with  positive  real  part  in  the  A  plane  for 
the  first  and  second  branches  of  even  and  odd  eigenvalues,  respectively,  without  including 
the  complex  conjugate.  Also  we  do  not  include  graphical  displays  for  the  axi-symmetric 
case  since  it  is  similar  to  the  odd  case.  These  graphics  have  been  obtained  for  some 
range  of  R  which  are  different  for  each  curve.  For  the  even  eigenvalues  the  first  curve, 
corresponding  to  the  complex  branch,  was  made  for  0  <  R  <  10°  while  the  second  branch 
for  0  <  R  <  100  .  For  odd  eigenvalues  the  first  curve  corresponds  to  0  <  R  <  90  .  and  the 
second  for  0  <  R  <  100  .  These  differences  arise  because  we  wanted  to  illustrate  the  cycles 
for  each  branch.  All  the  computations  were  done  on  the  VAX  11/780  at  the  Mathematics 
Research  Center  at  the  University  of  Wisconsin-Madison. 

Similar  to  4  we  present  several  tables  which  include  the  computation  involving  the 
eigenvalues.  For  R  =  0  ,  in  the  channel  case,  the  eigenvalues  may  be  computed  by  solving 
the  transcendental  equations  (2.8a)  and  (2.86)  .  Results  concerning  the  solutions  of  these 
equations  can  be  found  in  [6,  and  '11  .  However,  for  the  axi-symmetric  case  since  there  are 
not  previous  results  we  may  start  by  computing  the  first  few  eigenvalues,  with  positive  real 
part,  which  correspond  to  solutions  of  equation  (3.12).  As  we  said  in  §3  to  solve  equation 


(3.12)  we  used  solutions  to  (3.15)  as  initial  iterates  for  the  Newton  method.  At  the  same 
time  to  solve  (3.15)  we  used  the  New'ton  method  and  as  initial  iterates  the  expressions, 


with  7n  =  ?n  =  1,2-*-,  since  these  expressions  are  asymptotic  solutions  to 

(3.15).  In  each  case  we  have  found  that  only  two  or  three  iterations  were  needed  for  the 
Newton  method  to  find  solutions  to  (3.12)  and  (3.15)  with  six  decimal  places  of  accuracy. 
The  first  several  solution  to  the  equation  (3.12)  are  given  in  Table  1.  Also  we  are  including 
only  those  tables  for  w'hich  we  have  obtained  results  that  are  not  reported  in  4]  (the 
complete  set  of  tables  regarding  all  the  numerical  results  can  be  found  in  1  ),  for  example 
in  Tables  2,  3,  4,  5,  6  we  show  real  and  complex  eigenvalues  for  the  axi-symmetric  case, 
in  Tables  8,  9  we  give  the  critical  Reynolds  numbers,  and  eigenvalues,  corresponding  to 
graphics  (1.1a)  and  (1.16),  in  Tables  10,  11  we  show  the  eigenvalues  for  high  Reynolds 
number  calculations. 


As  pointed  out  in  A),  for  large  values  of  R  the  eigenvalues  which  have  large  modulus 
tend  to  be  less  accurate.  Since  the  complex  eigenvalues  with  negative  real  part  have 
larger  modulus  than  those  with  positive  real  part  they  may  be  a  sensitive  quantity  in 
the  computation.  We  computed  the  first  branch  of  odd  complex  eigenvalues  for  R  >  10 
with  different  tolerances  (the  tolerance  tol  is  a  parameter  in  PASVA3  that  controls  the 
grow  of  the  estimate  errors).  We  have  found  that  the  accuracy  of  the  eigenvalues  increased 
considerably  while  we  decreased  the  tolerance,  e.g.,  for  tol  =  10“°  and  10“ p  the  eigenvalues 
agree  up  to  six  decimal  places  for  Reynolds  numbers  in  the  range  10  <  R  <  50  and  they 
have  at  least  three  decimals  places  in  common  for  R  >  50  .  while  for  tol  =  10~8  and 
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10  10  the  eigenvalues  agree  up  to  six  decimal  places  for  Reynolds  numbers  in  the  range 
10  <  R  <  10000  .  We  present  the  corresponding  results  in  Table  7. 

6.  Eigenvalues  at  Critical  Reynolds  Numbers 

As  we  pointed  out  before  this  method  fails  to  compute  eigenvalues  at  critical  Reynolds 
numbers.  In  [4’  Bramley  and  Dennis  computed  some  eigenvalues,  that  they  suggested 
occur  at  critical  Reynolds  numbers,  e.g.  at  R  —  6.3  they  computed  the  eigenvalue  A  6.0, 
which  corresponds  to  the  first  critical  eigenvalues  of  the  second  odd  branch  (see  figure 
1.1b).  We  have  found,  according  to  our  computation,  that  R  =  6.3  is  not  exactly  a  critical 
Reynolds  number  but  it  is  very  close  to  it,  i.e.,  we  computed  two  eigenvalues  for  that  R,  one 
corresponding  to  the  increasing  solution  with  A  =  6.01267  (which,  apparently,  agrees  with 
the  value  given  by  Bramley  and  Dennis)  and  the  other  corresponding  to  the  decreasing 
solution  i.e.  A  =  5.85346.  A  similar  situation  happens  for  a  critical  eigenvalue  reported 
in  7  by  Gillis  and  Brandt.  They  computed  the  eigenvalue  A  =s  2.632,  presumably  a 
critical  one,  at  R  -  8.461  .  We  compute  eigenvalues  A  =  2.62085  and  A  =  2.63875  at  that 
Reynolds  number.  We  believe  that  it  is  difficult  to  determine  exactly  a  critical  Reynolds 
number  and  its  critical  eigenvalue.  Here  we  introduce  a  device,  which  may  be  useful,  to 
compute  better  approximations  to  critical  Reynolds  numbers  and  their  eigenvalues.  First 
we  assume  that  the  eigenvalues,  near  a  critical  Reynolds  number,  may  be  represented  by 
the  following  formula, 

A  =  Ac  a\  R  -  Rc  ,  (6.1) 

where  Rc  is  the  critical  Reynolds  number,  Ac  the  critical  eigenvalue  and  a  is  a  constant. 
For  Reynolds  numbers  /?i ,  /?2.  R3,  near  Rc  we  can  rewrite  (6.1).  after  some  algebraic 
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manipulation,  as 

Ri  =  R\  +  j4(A,  -  Aj)  +  B( A,  -  Aj)2  ,  (6.2) 

with  .4  and  B  given  by  A  =  ( 1  /a) 2  and  B  =  (l/d)2(2Aj  -  2AC)  respectively,  and  i  =  2,3. 
Then,  using  the  values  of  Rt,  A,  ,  we  solved  the  system  (6.2)  in  terms  of  A.B  .  We  can 
compute  Ac  by  the  formula 

Ac  —  A  i  -  —  ,  (6.3) 

while  Rc  is  obtained  using  (6.1).  This  approach  seems  to  work  well,  whenever  the  eigen¬ 
values  are  close  to  Ac.  In  Tables  8  and  9  we  show  all  the  critical  eigenvalues,  corresponding 
to  Figures  1.1a  and  1.1b,  computed  by  this  approach. 

7.  Asymptotic  Behavior  for  Large  R 

It  is  interesting  to  see  the  behavior  of  these  eigenvalues  for  high  Reynolds  numbers.  By 
high  Reynolds  number  we  mean  that  the  system  (4.1)  becomes  a  stiff  system  of  ordinary 
differential  equations.  In  [12],  Wilson  developed  a  theory  based  on  a  singular  perturbation 
analysis  to  investigate  these  eigenvalues.  He  considered  equation  (2.5),  corresponding  to 
a  channel  flow  problem  in  the  stream  function  formulation.  For  high  Reynolds  numbers 
Wilson  found  that  the  eigenvalues  approaching  zero  downstream  fell  in  two  categories, 
one  in  which  the  eigenvalues  approached  zero  as  0(/?_1)  and  in  the  other  as  0{R~ A)  . 
From  Tables  5  we  see  that  the  eigenvalues  behaving  like  0(R_1)  correspond  to  positive 
real  decreasing  eigenvalues,  which  according  to  Wilson,  these  eigenvalues  may  be  written 
as  A  —  \0R~l  -f-  0(R~Z)  with  A0  being  a  constant.  By  neglecting  smaller  terms  we  have 
that  A0  =  A R  .  Then  it  is  possible  to  use  the  numerical  data  to  compute  A0.  We  did  that 


for  R  =  105  and  the  Ao’s  that  we  found  agree  with  those  given  by  Wilson,  (see  [12],  Table 
3). 

The  second  type  of  eigenvalues  approaching  zero  downstream,  according  to  our  compu¬ 
tation,  are  the  even  complex  eigenvalues  (first  branch).  By  assuming  that  these  eigenvalues 
can  be  written  as 

A  =  XoR~^  i -  ,  (7.1) 

we  have  that  for  R  large  enough,  neglecting  smaller  perturbations,  we  may  compute  the 
Ao’s  ,  using  the  numerical  data,  by  Aq  «  A R^  .  In  Table  10  we  present  the  corresponding 
results. 

Notice  that  Ao  tends  to  a  constant  as  R  increases,  it  seems  to  be  that  R  has  to  be 
quite  large  before  neglecting  any  smaller  perturbation  in  the  asymptotic  expression  for  A. 
Unfortunately,  Wilson  did  not  mention  how  to  obtain  the  smaller  order  terms  in  (7.1). 
Using  the  numerical  data  we  may  attempt  to  compute  the  second  lower  order  term  in  (7.1) 
by  assuming  that  A  can  be  expanded  as 

A  =  bR~^  -f  cR~d - ,  (7.2) 

with  6.  c,d  constants  to  be  determined.  We  have  then  that  for  four  consecutive  eigenvalues 
Ai,A2,A3  and  A4  .  at  Reynolds  numbers  of  the  form  Rk  =  1 0/2^ _ x  .  the  following  formuia 

10~'a2  -  a3  =  Re(c)10'(i(10'i  -  ]0~d)R;d  .  (7.3a) 

and 

10~  ^  ct3  —  q4  =  /?e(c)  10 ~ 2<^(  1 0 -  '  -10 ~d)R;d  .  (7.36) 
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where  A*  =  a*  +  i/3*,  i.e.,  we  have  that  d  is  given  by 

d  =  log  f  10  _  .  °2 — — )  /log  10,  (7.4) 

MO  ■  q3  -  a4'/ 

while  b  can  be  determined,  using  again  (7.2),  by 

.  10-dA3-A4 

®  ,  _  1  * 

{I0~d  -  10"  t  )Rl  ■ 

and  c  is  determined  using  (7.2)  with  A  replaced  by  A4  . 

We  found  the  following  values  for  6,  c  and  d,  b  =  1.6400—0. 7813i,  c  =  —  0.1655-0. 0435i 

and  d  =  0.34676.  We  used  A*,  corresponding  to  R =  103,  104,  105,  10°.  Therefore  by 

substituting  these  values  in  (7.2)  it  is  possible  to  obtain  a  better  asymptotic  prediction. 

Notice  that  d  is  approximately  i.e.,  the  eigenvalues  behave,  approximately,  as 

A  =  (1.64  4-  0.7813i')/?_i  ^  0{R~“)  .  (7.5) 

In  Table  11  we  show  a  comparation  between  the  numerical  results  of  the  eigenvalue 
calculations  and  the  asymptotic  results  using  formula  (7.2).  The  agreement  is  satisfactory, 
notice  that  for  the  Reynolds  numbers  used  to  derive  formula  (7.2)  the  agreements  are 
better. 

8.  Conclusion 

The  eigenvalues  governing  the  rate  of  decay  of  perturbations  of  Poiseuille  flow  have 
been  studied.  The  procedure  to  compute  the  eigenvalues  is  based  upon  expressing  the 
eigenvalue  system  as  a  two-point  boundary  value  problem.  The  eigenvalues  are  computed 
one  at  a  time  jointly  with  the  eigenfunctions,  the  results  presented  are  in  agreement  with 
previous  computations  [4  and  12  .  For  high  Reynolds  numbers  we  show  that  the  first 


(7.5) 
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branch  of  even  eigenvalues  has  an  asymptotic  representation  of  the  form  (7.5).  We  think 
that  our  procedure  is  simpler,  more  direct  and  efficient  than  previous  methods  and  can 
be  used  with  relative  ease  on  a  computer.  Due  to  the  high  accuracy  obtained  on  using 
this  method  we  strongly  recommend  it  for  computations  involving  eigenvalue  systems  of 
ordinary  differential  equations. 
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complex  eigenvalues  with  positive  real  part  for  0  <  R  <  100 


4.46620- 1.467471  7.69410-1 .72697i  10.87457-1.89494 

4.29200— 1.42646i  7.52551~1.69447i  10.70765-1.86725 

3.76255-1.06197i  6.95388-1.31088i  10.11519-1.47824 

4. 24793-0. 78668i  7.31772-1.109201  10.40317-1.28557 

3.70720— 0.50937i  7.77552-0.31462:  10.82173-0.95050 

3. 79197-0. 45240i  7. 44769-0. 69724i  10.31659-0.90353 

4.17041— 0.27261i  7.48688-0. 33362i  10.22035-0.64556 


TABLE  3 


Complex  eigenvalues  with  smallest  positive 'real  part  for  R  >  100 


R 

— 

Ae(even) 

Ap(odd) 

A(axi-symm) 

250 

2.78306— 0.27633i 

0.72120— 0.37930i 

3.80091 -0.06650i 

500 

2. 71137-0. 15264i 

0.65637-0. 33642i 

3.97900-0. 17319i 

1000 

2.70753—0. 15583i 

0.59650— 0.29956i 

3. 86180-0. 15698i 

2500 

2.67659-0. 12999i 

0.38104-0.18322i 

3.50706— 0.09679i 

5000 

2.59490— 0.08062i 

0.34571— 0.16572i 

TABLE  4 

Axi-symmetric  complex  eigenvalues  with  negative  real  part  for  0  <  R  <  10 


!  R 

A 

A 

A 

0.25 

-4.5126— 1.4738i 

-7. 7378-1. 7309i 

-10. 9174-1.  S976i 

0.50 

-4.5601  — 1.4784i 

-7. 7821-1. 7331i 

-10.9607- 1.8986i 

1.00 

-4.6584  — 1.4823i 

-7.8727-1. 7318i 

-11.0486-1. 5949i 

2.50 

-4.9823— 1.4417i 

-8.1593— 1.6744i 

-11.3231-1. 8331i 

5.00 

-5.6481  — 1.0772i 

-8.6844- 1.2420i 

-11.8148-1. 4338i 

o 

o 

o 

-8.3161  — 1.5370i 

-1 1 .45  i  5—  1 .42  i  <i 

-14. 4887-1. 5590i 

2.84882 

1.23820 

0.63541 

0.32000 

0.12827 

0.06415 

0.03208 

0.00642 


1.86985 

0.95032 

0.38230 

0.19134 

0.09568 

0.01914 


1.89301 

0.76379 

0.38250 

0.19133 

0.03827 


3.21182 

1.27228 

0.63760 

0.31900 

0.06381 


1.90837 

0.95652 

0.47867 

0.09576 
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TABLE  7 


Odd  complex  eigenvalues  with  negative  real  part  for  R  >  10 


10 

25 

50 

100 

250 

500 

1000 

2500 

5000 

10000 


A0  (tol=10~G) 

-10.458682-  1.694815i 
-16.941646-2. 76307li 
-24.288883  — 4. 070337i 
-34.664217  — 5. 9667S9i 
-55.230448— 9. 7 108 19i 
-113.589146-11. 159176i 
-111.  168228-19. 879302i 
-176.178035  — 3 1.69208  li 
-362.979226— 36. 29621 3i 
-352.923824 -63. 768959i 


A0  (tol>  10“*) 

-10.458682-1.6948151 
-16.941646-2. 763071i 
-24.288883  — 4. 070337i 
-34.664089— 5. 966881  i 
-55.230451— 9.7 10823i 
-1 13.589146—  1 1 .1591 76i 
-11 1.168232— 19. 879303i 
-176.178035— 3 1.69208  li 
-362.97921  l-36.296142i 
-352.923830-63. 768960i 


TABLE  8 

Critical  Reynolds  numbers  and  critical  odd  eigenvalues 


««‘.o 

A? 

1  Ac 

I  8.4606 

2.63126 

6.2983 

5.89598 

9.1576 

3.42632 

6.8691 

!  6.97998 

25.6282 

2.67267 

17.1490 

5.91199 

!  26.3405 

3.15830 

18.0316 

6.76558 

50.8847 

2.70372 

31.6603 

5.95620 

51.3489 

2.95271 

32.6501 

6.63307 

84.1917 

2.74161 

49.3053 

6.00915 

50.7706 

6.51354 

71.4975 

6.19109 

TABLE  9 


Critical  Reynolds  numbers  and  critical  even  eigenvalues 


Ko 

K 

6.8846 

4.23405 

7.5805 

5.25272 

19.5660 

4.29246 

20.5465 

5.04126 

37.1756 

4.36997 

38.2524 

4.87873 

59.8274 

4.37883 

60.7830 

4.81713 

87.3756 

4.38935 

TABLE  10 


Even  complex  eigenvalues  for  high  Reynolds  numbers 


Ao 


I  1000  j  0. 59650—0. 29956i  1.60023-0.80362i 

I  10000  ;  0.43318— 0.20964i  1.61472-0. 78146i 

!  100000  j  0.31358-0. 15003i  1.62420-0. 77709i 

I  1000000  0.22651—0.108191  1.63013-0. 77865i 
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TABLE  11 


Numerical  vs  Asymptotic  Eigenvalues  at  high  Reynolds  numbers 


* 

Ae  (numer) 

Ae  (asympt) 

1000 

0. 59650-0. 29956i 

0.59625  — 0.28725i 

10000 

0. 43317-0. 20964i 

0.43317— 0.20780i 

50000 

0.34571  — 0. 16572i 

0.34572— 0.16552i 

100000 

0.31358-0.15003] 

0.31358— 0.15004i 

500000 

0.24986-0.  H934i 

0.24986-0.1 1950i 

750000 

0.23592-0.  Il269i 

0.23592-0.1 1271  i 

1000000 

0.22651-0. 10819i 

0. 22651-0. 10819i 
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